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The coherence properties of degenerate Bose gases have usually been expressed in terms of spatial correlation 
functions, neglecting the rich information encoded in their temporal behavior In this paper we show, using a 
Hamiltonian classical-field formalism, that temporal correlations can be used to characterize familiar properties 
of a finite-temperature degenerate Bose gas. The temporal coherence of a Bose-Einstein condensate is limited 
only by the slow diffusion of its phase, and thus the presence of a condensate is indicated by a sharp feature 
in the temporal power spectrum of the field. We show that the condensate mode can be obtained by averaging 
the field for a short time in an appropriate phase-rotating frame, and that for a wide range of temperatures, 
the condensate obtained in this approach agrees well with that defined by the Penrose-Onsager criterion based 
on one-body (spatial) correlations. For time periods long compared to the phase diffusion time, the field will 
average to zero, as we would expect from the overall U(l) symmetry of the Hamiltonian. We identify the 
emergence of the first moment on short time scales with the concept of U(l) symmetry breaking that is central 
to traditional mean-field theories of Bose condensation. We demonstrate that the short-time averaging procedure 
constitutes a general analog of the 'anomalous' averaging operation of symmetry-broken theories by calculating 
the anomalous thermal density of the field, which we find to have form and temperature dependence consistent 
with the results of mean-field theories. 

PACS numbers: 03.75.Hh 



I. INTRODUCTION 



The precise experimental characterization of the properties 
of Bose-condensed gases has motivated the development of 
theoretical methodologies that can provide accurate and com- 
prehensive descriptions of the condensed gas behavior. The 
fundamental theoretical framework is provided by many-body 
quantum field theory, but in general this becomes tractable 
only within approximation schemes, of which the most com- 
mon are based around Bogoliubov's idea of representing the 
condensed atoms by a classical mean field. In the very sim- 
plest form, this gives rise to the ubiquitous Gross-Pitaevskii 
equation, where the mean field is interpreted as the wavefunc- 
tion of the condensate. The solution of the Gross-Pitaevskii 
equation has provided a useful first approximation to a wide 
range of equilibrium and dynamical phenomena, but the equa- 
tion describes only the condensate, and neglects all sponta- 
neous and incoherent processes. There are many situations 
where the condensate is accompanied by a component of ther- 
mal atoms which can have an important influence on the sys- 
tem properties and behavior, and the early mean-field treat- 
ments have been extended to give some level of description 
of the noncondensed atoms, by employing factorization ap- 
proximations to the thermal component of the quantum field 
(U-O. Such self-consistent mean-field theories are built on the 
fictional (4) but convenient and intuitively appealing assump- 
tion that Bose condensation breaks the U(l) phase symmetry 
of the underlying quantum field Hamiltonian, resulting in the 
appearance of anomalous moments of the field: moments of 
the field operator [such as the mean field 0F(x))] which are 
formally zero in a state of fixed particle number, but which 



acquire nonzero values in the symmetry-breaking approxima- 
tion. These treatments have provided an improved descrip- 
tion of a range of equilibrium or near-equilibrium phenom- 
ena, but suffer from internal consistency problems, and have 
had limited success in describing the dynamics of the con- 
densate at higher temperatures (see Ref. [5 ] and references 
therein). We note that many of the equilibrium predictions of 
the symmetry-breaking mean-field descriptions are regained 
in more careful, number-conserving approaches (6^9), how- 
ever those methods have not provided a broadly tractable ap- 
proach for dynamical or higher temperature systems. 

In recent years, a set of techniques has been developed that 
provides a unified nonperturbative description of both equilib- 
rium and dynamical behavior of Bose gases for a temperature 
range from zero to close to the critical temperature. These 
so-called classical-field (or c-field) techniques fT0lfT2l have 
been used to provide a quantitative description of a number 
of key experimental results and regimes beyond mean-field 
theory (see ifTTTl for a summary of the broad range of recent 
applications). While the treatment superficially resembles the 
zero-temperature Gross-Pitaevskii theory, the interpretation of 
the central object of the theory, the classical field i/r(x), is 
very different: rather than the condensate wavefunction, it is 
an approximation to the Bose field operator, and provides a 
means to evaluate quantum-mechanical correlation functions 
and their time development. These correlation functions can 
be calculated by ensemble methods [Q21IT4), or for the case 
of equilibrium Bose-gas thermodynamics, by ergodic Hamil- 
tonian methods iTTSll . 

In this paper, we demonstrate that rich information is en- 
coded in the temporal behavior of Hamiltonian classical-field 
trajectories. Indeed in [16] we found that the temporal cor- 
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relations of a classical field revealed a strong signature of a 
quasicondensate-like structure in a spatially disordered (vor- 
tex liquid) phase. Here we consider the temporal correlations 
of a classical field containing a true condensate. The phase of 
a condensate is by definition 1 17 ] rigid across the spatial extent 
of the condensed mode, and the only condensate-phase fluc- 
tuations are global ones, which imply a diffusion of the phase 
over time (see Ref. fT8l and references therein). This diffu- 
sion restores the U(l) phase symmetry of the system in the 
ergodic (microcanonical) density of the field, and anomalous 
moments such as (^(x)> evaluated in this density therefore 
have vanishing values, consistent with the formal many-body 
theory for conserved particles. However, the time scale of this 
phase diffusion is typically long compared with the correla- 
tion times of thermal modes in the field fT8lL and we thus find 
that the condensate can be identified from the short-time aver- 
age of the field in a frame phase-rotating uniformly at the un- 
derlying (mean) phase-rotation frequency of the condensate. 
In this way phase- symmetry breaking emerges naturally from 
the Hamiltonian classical-field formalism. We demonstrate 
that this averaging procedure constitutes a general analog of 
the 'anomalous' averaging operation of symmetry -broken the- 
ories 1 1 1 by calculating the anomalous thermal density, which 
characterizes pairing correlations in the noncondensed com- 
ponent of the field which are induced by the interacting con- 
densate. 

This paper is organized as follows: In Sec. [TT] we give a 
brief outline of the equilibrium classical-field formalism we 
use here, and review its usual interpretation as a microcanoni- 



cal formalism. In Sec. Ill we discuss the emergence of a mean 
(first moment) of the classical field, and make a quantitative 
comparison to the condensate defined by the Penrose-Onsager 
measure of one-body coherence. In Sec. IV we consider 



the anomalous second moments which comprise the classical- 
field pair matrix, and construct the anomalous thermal density 
of the field. In Sec. [V] we summarize and present our conclu- 
sions. 



II. FORMALISM 

A. PGPE formalism 

The general formalism of (projected) classical-field meth- 
ods has recently been reviewed at length in ifTTl . but for 
the reader's convenience we will outline the projected Gross- 
Pitaevskii equation formalism we use in this work. The dy- 
namics we study are governed by the well-known classical- 
field Hamiltonian defined 



H cv = Jjx^(x)f7/ sp + ^| l A(x)| 2 ]^( X ), 



where the single-particle Hamiltonian is 

-h 2 v 2 



2m 



(i) 



(2) 



and the interaction strength £/ = 4nfi 2 a/m with m the atomic 
mass and a the s-wave scattering length. The projected clas- 
sical field is given by if/(x) = 2«gl a nY n (x), where the sum is 
over the finite set of single-particle eigenmodes [H sp Yk(x) = 
EkYkix)] with eigenvalues e n < E R , where E R is the single- 
particle cutoff energy. Defining the projector 

rm = Yj y «v> f d y Y *n(y)f(y)> e) 

we can express the Hamilton's equation for if/ix) obtained 
from Eq. ([I]) as 



= P ^ Hsp + W(x)| 2 )<A(x) 



(4) 



which is the projected Gross-Pitaevskii equation ifTTTl . The 
Hamiltonian H C f is invariant under the U(l) (gauge) transfor- 
mation xf/ix) — > i//(x)e l6 and has no explicit time dependence, 
so that the evolution described by Eq. ^ conserves both the 
normalization N c [i//] = J dx\i//(x)\ 2 of the classical field and 
the classical-field energy defined by H C f- In the microcanon- 
ical approach of the PGPE we follow here, finite-temperature 
equilibrium configurations of the classical field are obtained 
by evolving in real time randomized initial configurations con- 
structed with a particular energy E[ifs] = //cf[<AL such that the 
field naturally approaches thermal equilibrium, due to the er- 
godic nature |[T9l of the classical-field system. 



B. System parameters 

In the remainder of this paper we will specify quanti- 
ties in the characteristic units of the radial trapping poten- 
tial, quoting frequencies in units of a> r , distances in units of 
ro = ^h/moj r , times in units of oj~ 1 , and energies in units 
of %(jj r . We consider a system with oj z = ^oj r (represent- 
ing a typical three-dimensional trap geometry), and interac- 
tion strength N c Uo = V2 x 5QQftco r I r^. The corresponding 
ground (Gross-Pitaevskii) eigenstate of the system has energy 
E « 9N c tiQj r , and we choose the cutoff E R = 3lHoj r . We 
form random initial states EU Q21 with energies in the range 
E g [9.5, 24.0]N c Hoj r , which we allow to equilibrate by evolv- 
ing them in real time for a period of 120cj~ 1 , and perform our 
analysis on their subsequent evolution. 



C. Microcanonical interpretation 

Here we briefly remind the reader of the microcanonical 
(ergodic) interpretation of the PGPE applied to equilibrium 
systems (TT][T5l[20|. The method exploits the (empirical) fact 
that the PGPE trajectories are ergodic, and thus provide a sam- 
pling of the microcanonical density 



const //cf[*A] = E 
HqfW ± E. 



(5) 
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The trajectories ifs(x,t) cover the density P[\j/\E\ densely, 
and so averages in the density P[if/; E] are increasingly well- 
approximated by time-averages along trajectories if/(x,t) of 
increasing length. The implications of this for PGPE simu- 
lations are two-fold: First, a theorem due to Rugh [ 21 ] shows 
that the temperature of a microcanonical system can be ex- 
pressed as an average in its microcanonical density, and thus 
calculated from a time average. Second, equilibrium corre- 
lation functions of the classical field can similarly be defined 
as averages in the density ([5]), and thus evaluated from time 
averages. 

A correlation function of particular interest for characteriz- 
ing condensation in the classical field is the covariance matrix 
defined 

p(x,x') = (^(x)ifr(x% 

= ^ntftoXji*), (6) 

j 

[where (• • • ) J1 denotes a microcanonical average, i.e. an av- 
erage in the ensemble with density given by Eq. (|5J], which 
forms the classical-field analog of the one-body density ma- 
trix. In the second line we have used the fact that p(x, x') is 
Hermitian to express it in a diagonalized form, where the co- 
efficients {ft;}, indexed in order of decreasing magnitude, are 
the occupations of the corresponding modes faf;(x)}. By anal- 
ogy to the criterion of Penrose and Onsager (PO) ifTTlL con- 
densation in the field is signaled by the most highly occupied 
mode Xo( x ) having an occupation no which is significantly 
larger than all other occupations n,-. This definition in terms of 
correlations in the microcanonical density is an unambiguous 
measure of condensation in the simple equilibrium regimes in 
which it is applicable. Generalizations of this procedure based 
on short-time fluctuation statistics have been applied to more 
general scenarios, involving (e.g.) broken rotational symme- 
tries and nonequilibrium fields (T5J |22 - 24 j . In the remainder 
of this paper we will simply refer to the classical-field covari- 
ance matrix Eq. ([6]) as the one-body density matrix, and to the 
identification of its most highly occupied mode as the con- 
densate as the PO approach to quantifying condensation in 
the field. 

It is important to note that the microcanonical density 
[Eq. ([5])] inherits the invariance under gauge transformations 
i[r(x) — > ifr(x)e w of the Hamiltonian Eq. Consequently 
only the averages of quantities which are invariant under 
such transformations are nonzero in the microcanonical den- 
sity, which correspond of course to averages of operators 
which conserve particle number in the corresponding second- 
quantized field theory [25 ]. Sinatra and Castin fT8l have 
shown in a homogeneous geometry (where the condensate 
mode is a priori the k = plane- wave state), that the classical- 
field condensate undergoes a slow phase diffusion. This diffu- 
sion ensures that the gauge symmetry is restored in the micro- 
canonical density. In this paper, we show that on time scales 
short compared to the characteristic time scale of phase diffu- 
sion, the condensate is resolvable as the mean of the field in an 
appropriate frequency- shifted frame. Furthermore, we show 
that higher anomalous moments can be similarly defined in 



terms of short-time averages of fluctuations about this mean 
field. 



III. TEMPORAL COHERENCE 

A. Temporal coherence and the emergence of a nonzero first 
moment 

1. Identification of the first moment 

We begin by quantifying the coherence of the time- 
dependent field i//(x, t) via its temporal power spectrum, eval- 
uated at different spatial locations x (T6J. We define the tem- 
poral power spectrum of the classical field if/ at position x, 
evaluated over a period of length T as 

#(x;Q) = |g r {<A(x,0}| 2 , (7) 

where 3 r {/(0l denotes the Fourier coefficient taken from 
some arbitrary time origin 

% T {f(t)}=^ £ f(t)e iat dt. (8) 

In fT6l we applied this procedure to a classical field in a dis- 
ordered vortex-liquid state, in which spatial order of the sys- 
tem was strongly suppressed, and found a narrow peak in the 
power spectrum. The appearance of such a peak is consistent 
with analytical results obtained by Graham [26 ] which suggest 
that (quasi-)long-range spatial order of the Bose field is ac- 
companied by (quasi-)long-range temporal correlations which 
decay in a functionally equivalent way. Here we calculate the 
power spectrum for a classical field with the trapping and in- 
teraction parameters of Sec. |nj and energy E = l2.0N c hoj r . 
Using the PO approach, we find that this field exhibits a (true) 
condensate, with condensate fraction f c = no/N c = 0.70. We 
choose a sampling period of 40aj~ l , and approximate the in- 
tegral in Eq. ([5]) by a discrete sum over 1000 equally spaced 
samples of the classical field. In practice, we calculate the 
power spectrum at points in the z = plane, and average 
it over the azimuthal angle in this plane to smooth out fluc- 
tuations. We thus obtain the averaged power spectrum as a 
function of the radius r, which we present in Fig. [ija)- The 
oscillation frequencies Q we measure in the time-dependent 
field correspond, of course, to energies e = h£l in the quan- 
tum mechanical system. 

For comparison, on the same figure we also plot the pro- 
file of the harmonic trapping potential V(r)/h in this plane 
[parabolic green (gray) line] and the cutoff energy E R /h (hor- 
izontal line), and the classical turning point (vertical line) of 
the low-energy region L defined by their intersection. The 
most prominent feature in this plot is the strong peak in the 
power spectrum centered onQ = llAa) r , which is a signa- 
ture of the long-lived temporal phase coherence in the clas- 
sical field. We identify the frequency Aq of this peak as the 
condensate frequency . The broad, lower intensity background 
spectrum represents the thermally occupied excitations in the 
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FIG. 1 : (Color online) (a) Power spectrum H(x; Q) of the classical 
field on the plane z = (azimuthally averaged). Green (gray) lines 
indicate the trapping potential, cutoff energy and corresponding clas- 
sical turning point of the trap, (b) Space-integrated power spectrum 
N(Q) of the field [see Eq. (TT] )] as a function of the phase-rotation 
frequency £1 Parameters of the classical field are given in the text. 



classical-field system. It is worth pointing out that in the cen- 
tral region of the trap (r < 4r ), the background spectrum 
is strongly distorted by the presence of the condensate, with 
positive frequency components extending to approximately 
Aq + Er/H, and negative frequency components appearing with 
energies extending down to approximately Aq - E R /H l27l . 
Similar behavior was observed in Ref. (T6), and represents 
the restructuring of the excitation spectrum of the trap by 
an interacting (quasi-)condensate, which distorts the single- 
particle excitations of the system into the familiar Bogoliubov 
particle-hole pairs I71 I281E91 

The temporal coherence we observe results from the quasi- 
uniform phase rotation of the condensate: the phase of the 
condensate exhibits a uniform rotation at frequency Aq, super- 
posed with a slow diffusion. The width of the power spec- 
trum peak here is thus determined by the rate of this global 



condensate-phase diffusion. On time scales short compared 
with the characteristic time scale of phase diffusion, the con- 
densate has an approximately constant phase in a frame co- 
rotating at frequency Aq, i.e., short- time averages in this frame 
yield a nonzero first moment {if/) of the classical field. A key 
observation of this paper is that time averages constructed in 
this way are analogous to the anomalous averages which arise 
in symmetry -broken descriptions of Bose condensation Q, 
where the appearance of nonzero values for expectations of 
non-gauge-invariant quantities (i.e., the breaking of the phase 
symmetry) signals the presence of condensation in the field. 
We thus consider the classical field frequency- shifted by Q 



<A(x, t\ Q) = e iQt i//(x, 0, 



(9) 



and consider time-averages of this quantity formed from the 
same set of samples used to construct the power spectrum in 
Fig.[TJa). We define the time-averaged field 



0(x;Q) = (ifr(xj;Q)) t 
( = g r {^(x,0}), 



(10) 



where (•••)? denotes a time average over a given period T 
(40a;;: 1 in this case). The time-averaged field cp(x; £1) is there- 
fore the component of the classical field which phase-rotates 
like e~ lQ,t , and its norm square quantifies the total (i.e., space- 
integrated) power contained in the field at frequency £2, i.e., 



N(Q) 



= Jdx |0(x;D)| 2 = J dxH(x;Q). (11) 



In Fig. [TJb) we plot this power as a function of the fre- 
quency £2, and note that it exhibits a prominent peak at 
Q = 11.38<x> r . We identify the frequency at which the norm 
square of the time-averaged field (equivalently, the space- 
integrated power of the field) is maximized as the conden- 
sate frequency Aq, and the corresponding time-averaged field 
(p(x;Ao) as the classical-field condensate or mean field l30l . 
A nonzero time-averaged field occurs because the condensate 
has a reasonably well-defined phase on short time periods. We 
identify this quasi-definite phase as an analog of the conden- 
sate phase which emerges in symmetry-broken descriptions of 
Bose-Einstein condensation; in this view point, the first mo- 
ment 0(x; Aq) is the analog of the condensate wavef unction 
(^(x)) in such mean-field theories of Bose condensation. For 
notational convenience, we introduce the norm square of the 
mean field No = N(Ao), and the normalized mean-field mode 
function (po(x) = cp(x;Ao)/ VM)- The norm square A^o corre- 
sponds to the population of the mean-field condensate mode, 
and we indeed find No/N c = 0.706, in close agreement with 
the PO value for the condensate fraction (f c = 0.70). To fur- 
ther compare this temporal-coherence method of identifying 
the condensate with the PO approach, we calculate the over- 
lap of 0o ( x ) with the eigenvector Xo( x ) obtained by the PO 
procedure. We find 1 — K0oUfo)l ~ 1-4 x 10" 4 , i.e., the conden- 
sate orbitals obtained by the two different procedures agree to 
a very high accuracy. 
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FIG. 2: (Color online) Norm square N of the time- averaged field 
(i.e., space-integrated power of the classical field measured at the 
condensate frequency) as a function of the sampling period T. The 
solid line shows the value obtained from a single contiguous sam- 
pling of the classical field over period T. Circles (with error bars) 
show the mean (and standard deviation) of estimates No measured 
on 10 individual 400o;~ 1 sub-periods of the 4000o;~ 1 time series. The 
dashed line shows the expected (ensemble averaged) power, extrapo- 
lated from a least- squares fit to the means of the short- time estimates. 
Parameters of the classical field are given in the text. 



sample) period of the classical-field trajectory into 10 con- 
secutive sub-periods of length 4006L)" 1 (each of 10 4 samples), 
and regard these sub-periods as an ensemble of 10 distinct 
trajectories. For each member of the ensemble we calculate 
the power N(Aq\ T) as a function of T < 400&J" 1 . We then 
average over these 10 ensemble members to obtain a mean 
power estimate for each sampling period T. The means and 
standard deviations of these measurements are indicated by 
circles with error bars in Fig. |2| and by performing a least- 
squares fit of the expected power (N(Aq; T)) at the condensate 
frequency [Eq. ( |A4| ) in Appendix|A) to these mean power esti- 



mates, we estimate the phase-diffusion coefficient y « 10~ 4 &; r 
l3T1l . The dashed line in Fig. [2] extrapolates the expected be- 
havior of (N(Ao\ T)) to later times. Given this decay of the 
peak power with T, a rigorous estimate of the condensate pop- 
ulation would in principle be obtained by forming estimates 
(N(Aq, T t )) for multiple sampling period lengths T t , and ex- 
trapolating the resulting trend back to T = to estimate the 
'true' condensate population. However, due to the weak linear 
decay of the power spectrum peak at short sampling periods, 
we can accurately estimate the condensate population as the 
magnitude of the dominant peak in the power spectrum ob- 
tained over a short sampling period, for all but the smallest 
condensate fractions (see Sec. |IIIB 2| ). 



B. Dependence of the first moment on the field energy 



2. Temporal coherence and sample length 

The results obtained for the mean field have an important 
dependence on the averaging time. As discussed by Sinatra 
and Castin, the condensate phase exhibits diffusive evolution 
with time in the classical microcanonical ensemble 1 18 ]. Con- 
sequently, we expect the power in the classical field measured 
at the condensate frequency to decay with time, exhibiting a 
power-law tail N(Aq\ T) ~ 2/yT at long times, as discussed in 
Appendix|A| We illustrate this issue using the same simulation 
(E = l2N c hoj r ) as in the previous section. Increasing the sam- 
pling period to T > 50oj~ 1 the condensate frequency is more 
accurately resolved as Aq = ll.39oj~ l . We assume this value 
as a best estimate for the condensate frequency, and calculate 
the power at this frequency as a function of the measurement 
period T, up to a maximum measurement period of 4000^^. 
In Fig. [2] we plot the power measured at frequency Aq (solid 
line), and find that it decays in a nonuniform way with increas- 
ing T. However, the (normalized) mean- field orbital (po(x) we 
obtain at frequency Aq satisfies 1 - K0oUfo)l ^ 10~ 4 for all 
averaging periods T we consider; i.e., although the measured 
occupation of the condensate decays with increasing averag- 
ing period due to the diffusion of the condensate phase, the 
mode shape we obtain is relatively unaffected. The nonuni- 
form decay of the mean-field orbital 's occupation we observe 
is to be expected for a single trajectory, whereas we expect 
the scaling N ~ 2/yT to emerge from an average over a large 
ensemble of similarly prepared classical-field trajectories (cf. 
(T8ll ). It is possible, however, to infer y from a single tra- 
jectory, as we now show. We divide the total 4000a;;: 1 (10 5 - 



In the ergodic classical-field (PGPE) method, equilibrium 
field configurations of different temperatures can be formed 
simply by varying the (conserved) energy of the random ini- 
tial field configuration |20ll . In this section we investigate the 
behavior of the first moment introduced in Sec. lIII Al as the en- 
ergy (and thus temperature) of the classical-field equilibrium 
is varied, and compare its mode shape <po(x) and occupation 
No with the Penrose-Onsager condensate orbital Xo( x ) and oc- 
cupation no, respectively. We further compare the condensate 
frequency Aq to the microcanonical chemical potential p of the 
field obtained using the methodology of l2Tll32ll33ll . 



1. Condensate fraction 

We consider here the norm square No of the first moment 
III A| for various values of the 



0(x; Aq) defined as in Sec. _ 

classical-field energy E[if/]. In Fig. [3ja) we present estimates 
No for a range of classical-field energies, and compare them 
with the condensate occupations calculated by the PO ap- 
proach. The corresponding classical-field temperatures, cal- 
culated using the Rugh methodology (21] [32l [33), are also 
included in the figure. In practice, we calculated the PO con- 
densate by constructing the one-body density matrix [Eq. ([6])] 
from 3000 equally- spaced samples of the classical-field taken 
from a period of 1200GJ" 1 of the field evolution. We then di- 
vided this period into 30 consecutive sub-periods of length 
40oj~ 1 which we sampled at a higher resolution (1000 sam- 
ples per sub-period), from which we obtained 30 separate es- 
timates of the classical-field first moment 0(x; Aq). In each 
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FIG. 3: (Color online) (a) Condensed fraction of the classical field, 
as determined by the Penrose-Onsager procedure (circles) and by 
the time-averaging procedure (dots with error bars). Plusses indi- 
cate the microcanonical temperature of the field, (b) Discrepancy 
1 _ K0oUro>l between the (unit-normalized) first moment 0o( x ) and 
the most highly occupied natural orbital ^ ( x ) of the one-body den- 
sity matrix. 



sub-period we obtain the mean field as the time-averaged field 
of maximal norm, and we obtain (generally) distinct estimates 
of No, 0o (x) an d Aq from each series. In this way we exploit 
the ergodic character of the classical field to emulate sam- 
pling from an ensemble of similarly prepared trajectories (see 
Sec. |III A2| ). The red (gray) data points in Fig.[3ja) and their 
error bars represent in each case the mean and standard devia- 
tion of the norm squares of the 30 estimates of the mean-field. 
We observe that these estimates agree very closely with the 
PO condensate fractions f c (blue circles) throughout the range 
of energies presented. 

We also compare the mean-field orbitals obtained from the 
time-averaging procedure with the condensate orbitals ob- 
tained from the PO approach (see Sec. |III A In Fig. [3jb) 
we plot the the quantity 1 - |<0o[^o)l (averaged over the 30 
estimates) as a measure of the discrepancy between the two 
orbitals. We observe that for the energies presented (E < 
20N c fioj r ) the mean discrepancy is < 10" 2 . At higher ener- 
gies (corresponding to condensate fractions f c < 0.1), our 
temporal-coherence approach to identifying the condensate 
begins to break down: the mean-field orbital 0o(x) fails to 



match the PO orbital Xo( x ) (i-e. K0oUfo)l < 0.9) in an increas- 
ing fraction of estimates as the condensate fraction f c — > 0, 
and so for clarity we have not presented estimates for these 



energies here. This point is discussed further in Sec. IIIB 2 



2. Condensate frequency 

We consider here the dependence of the condensate fre- 
quency /io on the energy of the classical field. By our anal- 
ogy between the first moment of the classical field and the 
condensate wavefunction in mean-field theories (Sec. |III A| ), 
we associate this condensate frequency with the condensate 
eigenvalue appearing in such theories, which is itself closely 
related to the thermodynamic chemical potential of the degen- 
erate Bose-gas system (8j. In Fig.Qa) we plot estimates of the 
condensate frequency (red crosses), together with the thermo- 
dynamic chemical potential \i (blue circles with connecting 
line) of the classical field obtained from the Rugh analysis. 
At the very highest energies, we present results only for en- 
semble members for which our first moment analysis and the 
PO approach agree (i.e. 0o and^o overlap to within 10%). 
We observe that the condensate frequencies Aq and the chem- 
ical potentials \i agree very well for energies E < 20N c fioj r . 
Above this energy, the condensate frequencies Ao are consis- 
tently greater than the chemical potentials. This is expected 
behavior, as at a fixed total number of system particles, the 
two quantities differ by a factor of order 1 /N con a, where Af con( j 
is the condensate occupation (8j[34l. Davis et al. l35l ar- 
gued that equipartition of energy in the classical-field model 
predicts the relationship 



Mo 



k B T 
No ' 



(12) 



In Fig.[4ja), we plot the quantity Aq-IcbT '/HNo (black plusses), 
where the temperature T is that obtained from the method of 
Rugh, and thus find that our results are in reasonable agree- 
ment with the prediction of Eq. ( fl"2| ). 

We now consider how the total power spectrum of the clas- 
sical field varies as a function of the field energy. In Fig.|4jb) 
we plot the power spectrum N(Q), averaged over the 30 indi- 
vidual 40oj~ 1 sampling periods, for field energies in the range 
E e [9.5, 24]N c fioj r . At the lowest energies the behavior of 
N(Q) is as in Fig.[TJb): the function exhibits a prominent peak 
which we identify with the condensate, and a broad back- 
ground we associate with thermal excitations. As the energy 
(and thus temperature) of the field is increased, the conden- 
sate peak decays, and the 'wing' of thermal excitations grows 
until it is of the same magnitude as the condensate peak, and 
at the highest temperatures only the thermal background re- 
mains. This explains why our approach to identifying the 
condensate begins to fail as the temperature approaches the 
phase transition: although a temporally coherent condensate 
may still be present in the field, it becomes increasingly likely 
that the peak power in any particular estimate of the power 
spectrum corresponds instead to thermally occupied modes, 
which eventually swamp the condensate completely. 
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FIG. 4: (Color online) (a) Condensate frequency obtained from the 
time- averaging procedure (red crosses), and microcanonical chem- 
ical potential of the classical field (blue circles), for classical field 
equilibria with different energies. Black plusses plot the RHS of 
Eq. |T2| . (b) Norm square of the first moment as a function of the 
phase-rotation frequency of the frame in which it is constructed (i.e. 
space-integrated power spectrum of the field), for classical-field sim- 
ulations with different energies. 



IV. PAIRING CORRELATIONS 

In the previous section, we have identified that the conden- 
sate present in the classical field can be well-characterized as 
the time-average of the appropriately frequency-shifted field. 
We now show that more general anomalous moments of the 
field can be obtained from time-averages in the same phase- 
rotating frame. In this approach, condensation in the classical 
field is thus accompanied by the appearance of anomalous mo- 
ments of all orders, in direct analogy to the emergence of gen- 
eral anomalous correlation functions in symmetry-breaking 
accounts of Bose-Einstein condensation fT1l36ll. 

In terms of the Fock-space decomposition ^(x) = 
YjidiYi(x), the emergence of a mean field in the second- 
quantized formalism is equivalent to the appearance of 
nonzero first moments {(#/)}. The next- simplest anomalous 
averages, the quadratic moments {{didj)} (and their conju- 
gates), arise due to the effect of interactions which 'mix' 
the single-particle creation and annihilation operators to form 
quasiparticle operators b ~ ua + v*d^ |25 ]. Consequently, the 
occupation of quasiparticle modes results in the appearance 
of nonzero moments of single-particle operators of the form 



(didj), which represents correlations between pairs of parti- 
cles. Like the mean field itself, these moments are formally 
zero in a state of fixed total particle number, although anal- 
ogous quantities can be defined in particle-conserving terms 
151 . Due to the appearance of these pairing correlations, in 
order to accurately characterize the weakly interacting Bose 
gas and its excitations, one must consider not only the one- 
body density matrix pjj = (d^di), but also the pair matrix 
Kij = (djdi) E3. 

In the remainder of this section, we will demonstrate the 
application of our temporal averaging procedure to the evalu- 
ation of quadratic anomalous moments of the classical field: 
by estimating the pair matrix k(x, x') = (i//(x)i//(x')) of the non- 
condensed component of the field, we calculate the anoma- 
lous density which characterizes pairing correlations in the 
thermal component of the field. We note that signatures 
of such pairing correlations have been observed previously 
in classical-field calculations 1231 . where anomalous values 
gf^ = (\ai\ 4 ) / (\ai\ 2 ) 2 > 2 were obtained for the second-order 
coherence functions of density-matrix eigenmodes. We note 
also that a temporal signature of the anomalous density has 
previously been observed l37l in homogeneous classical-field 
simulations, in which the anomalous density is uniform. 



A. Methodology 

We seek here to characterize pairing correlations in the ther- 
mal component of the classical field, i.e., the component of 
the field orthogonal to the condensate |3], which is obtained 
by projecting out the condensed component of if/(x, t), i.e. 

ifs ± (x,t) = ifr(x,t)-(f)o(x) J dx f f (x f )i//(x\t)dx\ (13) 

It is important to note that we form ^(x, t) on a given (40cj~ 1 ) 
time period by projecting out the mean field obtained over 
the same period, so that (anomalous) averages constructed 
from ^(x, t) over this period are formed on the same foot- 
ing as the mean field itself. We transform ^(x, t) to the same 
phase-rotating frame as the condensate, forming \j/^(xj) = 
e^^ix, t), and then calculate the pair matrix 



K ± (x,x') = (^(x)i(f ± (x')) t . 



(14) 



The most well-known characterization of the anomalous cor- 
relations described by the pair matrix is given by the anoma- 
lous density [ 1 ] which we identify as the diagonal part of the 
pair matrix 



m(x) = #(#(x)) f E^(x,x). 



(15) 



We find that the general form of mix) is apparent from a single 
estimate of /c x (x, x'), over a temporal period AQco~ l . However, 
large fluctuations are present in such a single estimate, which 
is to be expected, as the correlations we seek to resolve here 
are rather subtle as compared, for example, to the coherence 
of the condensate. In order to resolve the anomalous density 
more clearly, we therefore average over multiple estimates of 
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m(x); i.e., for each of 30 consecutive 40qj~ 1 periods, we form 
both the mean field <po(x) and the corresponding anomalous 
density mix). The phase of mix) is only meaningful in rela- 
tion to the phase of the mean field <po(x) itself, and so for con- 
venience, we choose the overall phase of the classical field in 
each sampling period such that (po(x) is maximally real. This 
choice of the phase of mix) relative to a real and positive con- 
densate wavefunction corresponds to the traditional choice in 
mean-field theories. Forming multiple estimates of m(x) in 
this way allows us to calculate both the mean and the variance 



(a) 



of this quantity, as we show in Sec. IV B 
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B. Anomalous density 




y [units of ro 



The mean anomalous density calculated by the procedure 
described in Sec. IV A is mostly real and negative (i.e., has 
phase opposite to that of the mean field), in agreement with 
the results of mean-field theory calculations (3| [34j [38), but 
exhibits some small complex- valued fluctuations due to the fi- 
nite ensemble size. Denoting averages over estimates by an 
overbar, we plot in Fig. [^a) the negative -Re{m(x)} of the 
mean anomalous density on the z = plane, and the local 
standard deviation in estimates Sm(x) = [|m(x)| 2 - |m(x)| 2 ] 1/2 
of the anomalous density on this plane, calculated for a simu- 
lation with E = 14.5N c fiCL> r (for which the condensate fraction 
f c = 0.50). The anomalous density has the spatial structure 
expected from mean-field calculations (3] [34] [38): it resides 
primarily in the region where the condensate exists, and its 
absolute value exhibits a shallow 'dip' in the center of the 
trap. The standard deviation 8mix) indicates that the great- 
est variance in density estimates occurs around the (circular) 
maximum of |m(x)|, while much less variation occurs in esti- 
mates of the density in the central dip. 

The anomalous density shown here exhibits a very high de- 
gree of rotational symmetry about the z axis, but in general the 
anomalous density we obtain is distorted (the central 'dip' in 
its absolute value becomes saddle-shaped along some random 
axis). We identify this as a result of persistent center-of-mass 
(dipole) excitations of the field (39), which are 'frozen in' dur- 
ing the thermalization of the field. More generally one might 
regard the classical field as having condensed into an excited 
center-of-mass mode, and consider the correlations of the field 
in a frame following this motion [40]. In Fig.[5jb) we plot the 
azithumally averaged anomalous density on the plane z = 0, 
together with the similarly averaged densities of the conden- 
sate [A^ol0o(x)| 2 ] and the orthogonal thermal component of the 
field [|^(x)| 2 ], for comparison. We observe that the magni- 
tude of the anomalous density in the center of the trap is an 
appreciable fraction of that of the (normal) thermal compo- 
nent of the field, in agreement with Refs. | 3 , 341138). 



C. Dependence on field energy 

Finally, we consider the dependence of the anomalous den- 
sity on the energy (or equivalently, the temperature) of the 
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FIG. 5: (Color online) Anomalous density m(x) of the field, for case 
E = \4.5N c fta) r . (a) Shape of -Re{m(x)} on a slice through plane 
z = (upper surface), and standard deviation Sm(x) of anomalous- 
density estimates on the same plane (lower surface), (b) Azimuthally 
averaged density of the condensate mode (as determined by the time- 
averaging), complementary (orthogonal) thermal component of the 
field, and anomalous density, in the plane z = 0. 



projected classical field. In mean-field theories, the anoma- 
lous density (after any renormalization (§1 [34] [41)) becomes 
small as the temperature of the system approaches zero (due to 
the weak occupation of quasiparticle modes in this limit), and 
also as it approaches the critical temperature (due to the quasi- 
particle modes becoming more single-particle-like as the con- 
densate is depleted). This behavior is often cited as a justifica- 
tion for the neglect of the anomalous density in self-consistent 
theories (the so-called Popov approximation O) in these two 
limits. In order to characterize more fully the temperature- 
dependent behavior of the anomalous density, we follow (3) 
and calculate its integrated value M = Jdxm(x). In Fig. 
we plot the real part of M (neglecting a small imaginary part 
that arises from incomplete convergence of the averaging - 
see Se c.|IVB| ) as a function of the classical-field energy. As 
in Sec. IIIB| we present for each energy the mean and stan- 
dard deviation (error bars) of only those estimates obtained 
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FIG. 6: (Color online) Integrated value of the anomalous density 
J dx m(x) as a function of field energy. 



from averaging periods which produced an accurate conden- 
sate mode (Sec. |III A| . We observe that the behavior of M(E) 
is consistent with the results of mean-field theories (3] [38), 
with its absolute value \M\ reaching its maximum at interme- 
diate energies (temperatures), and rapidly decreasing as we 
approach both the zero-temperature and critical regimes. 

We note that the well-known issues of ultraviolet diver- 
gence of the anomalous density in mean-field theories arise 
from the zero-point occupation of quasiparticle modes (H EU 
|4H . which is of course not present in the classical-field model. 
Also, although the results of classical-field calculations are 
in general dependent on the cutoff energy, the contribution to 
the anomalous density from successive quasiparticle modes 
rapidly decreases with increasing energy of the modes, as 
the modes return to a single-particle structure. The require- 
ment (for our treatment) that the anomalous density is well- 
contained in the low-energy region (condensate band fiTTl ) de- 
scribed by the PGPE is thus precisely the requirement that the 
cutoff is effected at such an energy that the interacting Hamil- 
tonian has become approximately diagonal in the single parti- 
cle basis {Yk(x)} (satisfied in practice for Er/ji > 3 [42, 43 ]). 

Finally, we note that the maximum (absolute) value of the 
integrated anomalous density occurs when f c « 0.5, but re- 
mind the reader that this refers only to the proportion of 
the below-cutoff field which is condensed. Although the en- 
tire anomalous density should be well-described by the low- 
energy Hamiltonian PGPE dynamics, one would have to in- 
clude the contribution of above-cutoff atoms to the normal 
thermal density of the field [11] in order to draw quantitative 
comparisons with (e.g.) the mean-field theory calculations of 



0. 



V. CONCLUSIONS 

We have demonstrated that in the Hamiltonian PGPE the- 
ory, classical- wave condensation is accompanied by long- 
range temporal coherence limited only by the slow diffusion 



of the condensate phase. This gives rise to the appearance of 
a nonzero first moment of the field, as defined by short-time 
averages in an appropriate phase-rotating frame. We identi- 
fied the emergence of this moment with the concept of U(l)- 
symmetry breaking that is central to self-consistent mean-field 
theories. We showed that the mean field obtained by short- 
time averaging agrees well with the condensate identified by 
the standard Penrose-Onsager approach, except for close to 
the critical regime associated with the transition to the nor- 
mal phase. The condensate eigenfrequency obtained by this 
temporal analysis exhibits the behavior predicted for the con- 
densate eigenvalue in the most sophisticated mean-field ap- 
proaches [8], i.e., it agrees closely with the thermodynamic 
chemical potential at low energies and diverges away from 
it in inverse proportion to the condensate occupation as the 
phase transition is approached. By calculating the pair matrix 
and anomalous density of the noncondensed component of 
the field, we demonstrated explicitly that time averages in the 
frame rotating at the condensate frequency allow the calcula- 
tion of more general anomalous moments. We observed the 
anomalous density to exhibit the expected behavior (3), with 
its magnitude reaching its maximum at intermediate temper- 
atures and decreasing as both the T = and critical regimes 
are approached. 
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Appendix A: Phase diffusion of the condensate 

1. Single phase-diffusive mode 

Let us first consider a single-mode model of the condensate, 
in which the amplitude ao(t) of the condensate mode (with 
condensate frequency Ao) exhibits phase diffusion. We assume 
that the mode does not exhibit any number fluctuations, which 
is precisely the condition that g^ 2) = (\ao\ 4 ) t /(\ao\ 2 ) 2 = 1> 
which is well satisfied away from the critical regime I23ll44l . 
We thus have a = \a \e m , and defining (p(t) = 0(t) - 0(0), 
we assume var {<£>(£)} = {(fit) 2 ) - {(fit)) 2 = 2yt flU, where y 
is the (phase) diffusion coefficient, and (• • • ) denotes an aver- 
age over realizations of the amplitude ao(t) (i.e., an ensemble 
average). This is precisely the behavior of the Kubo oscilla- 
tor l45l stochastic process, which obeys the (Ito) stochastic 
differential equation 

dao = [(-iAo - j)dt + i j2ydW(t)]ao(t), (Al) 
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where dW(t) is a real Wiener increment, which satisfies 
(dW(t)dW(f)) = S(t - f)dt. By studying this simple model 
we hope to gain insight into the behavior of our diffusive con- 
densate mode. 

We consider the power spectrum of the mode obtained over 
a period T, 



i 1 r T 

Ar (0) (Q;D= - dt 
11 Jo 



e iQt a (t)\ 



(A2) 



This power spectrum is itself a stochastic process (developing 
in T), i.e., it varies between realizations of the oscillator. We 
therefore consider its mean (Af (0) (Q; T)). Using the known 
result (ao(t)aQ(s)) = \ao\ 2 exp[-iAo(t - s) - y\t - s\] B31L we 
find 



<Af (0) (Q; T)) 



1 



r r( r 2 + A 2 ) 



(A3) 



T 2 (y 2 + A 2 ) 2 
+[e~ yT cos(A T) - l](y 2 - A 2 ) - 2yA sin(A 



which we have written in terms of A = £1 - Aq for compact- 
ness. In the limit of no diffusion (y — > 0) we regain the re- 
sult N (0) (Q;T) = |a | 2 sinc 2 [|(Q - A )t] appropriate to the 
resolution of a single frequency by a measurement of finite 
duration T. In the limit of a measurement made on a time 
scale long compared with the characteristic diffusion time, i.e. 
yT » 1, we regain the Lorentzian spectrum of the Kubo oscil - 
latorAf (0) (ft;r) = (2\a \ 2 y IT)l[y 2 + (Q-A ) 2 ]. FromEq. ( [A3] ), 
the power measured at the underlying frequency Aq of the os- 
cillator can be obtained by setting A = 0, giving 



only by the time scales on which the power decays. Qualita- 
tive differences between the two types of mode thus only ap- 
pear in second-order (and higher) correlation functions, which 
are sensitive to number fluctuations. 



2. Multimode description 

In general the condensate mode is only one mode in a mul- 
timode field which contains other, thermally occupied modes. 
We expect the thermal field to be well described in the basis 
of Bogoliubov modes {(M,-,Vj)} orthogonal to the condensate 
mode [49], and thus assume 

<A(x, t) = a (t)xo(x) + Yj (bj(t)uj(x) + Z?*(0v*(x)), (A5) 

j 

where, to gain simple insight into our measurements of 
the field, we assume that the {bj(t)} are complex Ornstein- 
Uhlenbeck processes which are uncorrelated with one another 
and with the condensate. The total power spectrum of the field 
is thus 



(N (0 \Ao;T)). 



mi 



(yT)- 



{yT-(l-e-y T )\ 



(A4) 



For short time periods r « 1/y (such as we consider in 
the main text), the norm square of the mean field decays like 
~ 1 -yT/3, while at long times it decays like ~ 2/yT. It is 
important to note that this same functional form would be ex- 
hibited by (e.g.) a complex Ornstein-Uhlenbeck process l45ll . 
which one might reasonably assume as a model for a ther- 
mally occupied mode I46ll48l in a classical-field approxima- 
tion: the 'bare' (i.e. infinite sampling time) power spectrum 
of such a mode is similarly Lorentzian, and so we expect the 
same behavior both for two-time correlations [\(a* (t)a(0))\ ~ 
e~ yt ] and for the measured power N(Q T), and the two cases 
(i.e. condensate and thermal mode) are thus distinguished 



(N(Cl;T)) 



<Af (0) (Q;r)> J dx\xo(x)\ 2 (A6) 
+ ^<A^' ) (Q;r)> J dx\ Uj (x)\ 2 + |v 7 -(x)| 2 , 



where 



(N (j \£l;T)) 



dte iat bj(t)\ 



(A7) 



behave similarly to {N^\£l\ T)), except that they are centered 
on the frequencies e B .j% of the Bogoliubov modes, and atten- 
uate much more rapidly with T (yj » yo)- There is therefore 
power in the field at a range of frequencies, however, on times 
T » l/yj we have (N(Q;T)) « Af (0) (Q; T) f dx\xo(x)\ 2 and, 
moreover, 

(I J dte iA ^(xj)) = 1 ~^ t 7T \a \(e im ) X o(x), (A8) 

where the appearance of the expectation (e l6 ^) of the initial 
complex phase emphasizes that the condensate phase varies 
randomly between ensemble members, 'breaking' the U(l) 
symmetry in any particular realization. 
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